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WEINAN E, JIANFENG LU, AND XU YANG 

Abstract. The behavior of interacting electrons in a perfect crystal under 
macroscopic external electric and magnetic fields is studied. Effective Maxwell 
equations for the macroscopic electric and magnetic fields are derived starting 
from time-dependent density functional theory. Effective permittivity and 
permeability coefficients are obtained. 



1. Introduction 

This paper is a continuation of our study on the macroscopic behavior of in- 
teracting electrons in a crystal. In the previous paper [8], we studied the Bloch 
dynamics of a single electron in a crystal and introduced the Bloch- Wigner trans- 
form for studying the semi-classical limit of Schrodinger equation. We also gave a 
simplified derivation of the Berry curvature term in the effective dynamics. In this 
paper, we study the collective behavior of the interacting electrons in an insulating 
crystal under applied electric and magnetic fields. We derive the effective Maxwell 
equations in this case using systematic asymptotics. In particular, we obtain the 
effective permittivity and permeability coefficients for these materials. 

From a macroscopic viewpoint, the behavior of crystals can be characterized as 
follows: 

(1) Mechanically, crystals respond to applied stress by deforming the crystal 
lattice. 

(2) Crystals respond to applied electric and magnetic fields by distorting the 
charge-spin distribution, or by motion of free electrons. This generates 
electro-magnetic responses. 

The mechanical and electro-magnetic responses can be coupled together, generating 
piezo-electric, magnetorestrictive and ferro-elastic effect, etc. The main purpose of 
this series of work is to provide a systematic understanding of these macroscopic 
phenomena and derivation of the effective macroscopic models from "first princi- 
ples". 
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As the first principle, we choose to work with the density functional theory 
p~3lfT4]fT9| instead of the many-body Schrodinger or Dirac equations. This is because 
that density functional theory has proven to be extremely successful for the kind 
of issues we are interested in, and is at the present time the only tractable and yet 
reliable models for electronic matter. Here by density functional theory, we mostly 
mean Kohn-Sham density functional theory that rely on orbitals, as is done in this 
paper. But occasionally we also resort to orbital- free density functional theory, 
such as the Thomas-Fermi type of models, to illustrate some of the issues. We refer 
to [Tll2ll4l[6l lT5Hl8] for the mathematical works done on density functional theory. 
Closely related are the works on Hartrcc or Hartree-Fock models, which have also 
been used as the starting point for analyzing the behavior of crystals. 

When the crystal is elastically deformed, continuum mechanics model can be 
derived from the Cauchy-Born rule (extended to electronic structures). This was 
done for the Thomas- Fermi- von Weiszacker model in [4] . In a series of works by 
E and Lu [T0lU2| . the Cauchy-Born rule was validated for nonlinear tightbinding 
models and Kohn-Sham density functional theory. One of the important ingredi- 
ents in these works is the identification of sharp stability criteria when the model 
has exchange-correlation energy which might be non-convex. The issue of stability 
does not occur in Thomas-Fermi- von Weiszacker, Hartree or reduced Hartree-Fock 
model, since these models do not include exchange-correlation energy. Further in 
this direction, E and Lu studied in [5] the continuum limit of the spin-polarized 
Thomas-Fermi-von Weiszacker-Dirac model under external macroscopic magnetic 
fields. Under stability conditions for plasmon and magmon, a micromagnetics en- 
ergy functional was derived. 

One interesting by-product of the work in is an effective model for the 
macroscopic electric potential as a result of the crystal deformation, which exhibits 
a coupling between the mechanical and electric responses. 

Cances and Lewin studied the reduced Hartree-Fock model for a crystal under a 
macroscopic external potential and proved that the implied macroscopic potential 
satisfies an effective Poisson equation. In particular, they established the validity 
of the well-known Adler- Wiser formula for the permittivity tensor [5] . 

In this work, we consider the time-dependent Kohn-Sham density functional 
theory in the presence of external macroscopic electric and magnetic fields. The 
questions of interest are whether macroscopic Maxwell equations that describe the 
electromagnetic fields can be derived in the continuum limit from the underlying 
microscopic theory, and in particular, how to obtain effective permittivity and per- 
meability for materials from electronic structure models. We resolve these issues 
using asymptotic analysis. To rigorously justify the asymptotic derivation, one 
needs to identify correct stability conditions for time dependent models. This will 
be left to future publications. 
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The paper is organized as follows. In Section^ we introduce the time-dependent 
Kohn-Sham density functional theory. Section [3] describes the model setup and 
presents the main results. The asymptotic derivation is given in Section 31 Section 
[5] and Section [6l We make conclusive remarks in Section [JJ 

2. Time-dependent density functional theory 

Time dependent density functional theory (TDDFT) [19] is an extension of 
(static) density functional theory to the dynamics of interacting electrons. In 
TDDFT, the electron dynamics is governed by N one-electron time dependent 
Schrodinger equations with effective one-body Hamiltonian depending on electron 
density and/or electron current densityQ 

The TDDFT model takes the following form in physical units in R 3 , 
dijj 1 



(2.1) ih^- = (-ihV -~ c (A + A cxt )) 2 ^ + e(V + V cxt )^, 



e 



(2.2) -A0=-(p-m), 

cot \cot J ce 

(2.4) V • A = 0, 

(2.5) V(t,x) = cf>(t,x) + ri{p(t,x)). 

Here ifij, j = 1, . . . , N, is the one-particle wave function, A is the vector potential 
and (j) is the scalar potential generated by electrons. The electric and magnetic 
fields are given by 

BA 

E=-V(f>-—, B = WxA. 
at 

The system is invariant under the gauge transform, 

A^A + Vx, 0^^-^, 

and hence we fix the Coulomb gauge (j2.4|) in the model. A ext and V ex t are the 
external vector and scalar potentials. The electron number density and electron 
current density are denoted by p and J respectively in the equations, and are given 
in terms of {ipj}f = i by 

N 

(2.6) p(t,x) = ^|^(i J K)| 2 , 

j=i 

h N e 

(2.7) J{t,x) = — V 1, 3m (ipJ(t,x)Vipj (t,x)) —p(t,x)A(t,x). 

m e ^— ' J m e c 



^When the effective Hamiltonian depends on electron current density, the model is usually 
called time dependent current density functional theory (TDCDFT) [201121] in physics literature, 
although we still use the name of time dependent density functional theory in this paper. 
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The function m(x) is the background charge density contributed by the nuclei. We 
assume that the nuclei are fixed so that m(x) is independent of time. In (|2.1j) . we 
have the physical constants electron mass m e , electron charge e, Planck constant 
h, dielectric constant in vacuum eo and speed of light in vacuum c. 

The electric and magnetic fields are given by vector and scalar potentials (in the 
Coulomb gauge), 



We make some remarks about the model. 

(1) The spin is ignored in the above TDDFT model. As a result, only the orbital 
magnetization is considered, while spin magnetization is not present. The 
extension to include spin in the model is straightforward though. 

(2) We adopt the adiabatic local density approximation [71122] for the exchange- 
correlation potential, denoted as rj in equation (|2.5I) . This means that the 
exchange-correlation potential is a function of local electron density only. 
No exchange-correlation vector potential is included in the model. Gen- 
erally, the exchange-correlation potential can depend on the local electron 
current density and the derivatives of electron density and current density. 
Exchange-correlation vector potential can also be added. The extension 
to these general models is in principle possible, but will complicate the 
formulations and derivations in the discussions below. 

(3) The model agrees with what physicists commonly use in practical applica- 
tions (for instance [3]). Of course, whether the model gives a good predic- 
tion of the time-evolution of electronic structure depends on the choice of 
pseudo-potential, the choice of exchange-correlation functional, and some- 
times requires additional terms like exchange-correlation vector potential. 
We will not go into the details of this discussion. 

Nondimensionalization and high frequency scaling 

We consider the situation when the applied external fields to the system have 
a much larger characteristic length compared to the atomistic length scale (lattice 
parameter). For this purpose, we perform nondimensionalization to the set of 
equations and identify small parameters. 

We introduce two sets of units to rescale the system. One is the microscopic unit 
in which we denote the units of time, length, mass and charge as [t], [I], [to], [e]; 
the other is the macroscopic unit in which we denote the units of time and length 
as [T], [L]. It means that for example the characteristic time scale for macro- 
scopic fields is [T], while that for microscopic fields is We will consider the 
macroscopic behavior of the system under macroscopic external potentials within 



(2.8) 



E 








(2-9) 



B = V x A. 
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the high frequency regime, in other words, the regime 

[£]»[/]. 

The small parameter is identified as e = Physically, the high frequency 

regime means that we are interested in the dynamics of the electronic structure and 
the corresponding dynamics of the electromagnetic fields on the time scale that is 
comparable to the characteristic time scale of the quantum system. At longer time 
scale, different physical phenomena might occur and is not covered by the results 
here. In particular, this is different from the scaling used when considering the 
semi-classical limit. 

Using these two sets of units, we can represent all physical constants and quanti- 
ties in suitable units so that they become nondimensional and have values of order 
0(1). For example, Planck constant, vacuum dielectric constant and speed of light 
can be written as 

* i v NW 2 , H 2 W 2 , v [L] 

h = lx —> £o = lx W c = lx [7y 

The temporal and spatial derivatives are rescaled as 

d_ J_d_ 1 

dt ~~ * [T]dt } ~^ [L] ' 

The physical quantities are rewritten as 

where A, V, p , J are nondimensional quantities. 

Substituting all the above into the system (12.1l) - (|2.5j) produces the nondimen- 
sionalized TDDFT equations (the tildes are dropped for simplicity), 

(2.10) i ^T = \ H £V - < A + ^ext)) 2 ^ + (V + V eKt )iP j} 

(2.11) - A(/> = e(p - m), 

(2.12) ^A-AA+^V^J, 

(2.13) V-A = 0, 

(2.14) V(t, x) = (p(t, x) + r,{e^p{t, x)). 

The density and current are given by 

N 

(2.15) p(t,x) = Y / \^(t>v)\ 2 i 

N 

(2.16) J(t,x) = e^2Jm(ip^(t,x)Wipj(t,x)) - sp(t,x)A(t,x). 
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3. The effective Maxwell equations in crystal 

3.1. Unperturbed system. Let L be a lattice with unit cell T. Denote the re- 
ciprocal lattice as L* and the reciprocal unit cell as T* . We consider system as 
a crystal eL, so that e is the lattice constant (the micro length scale used in the 
non-dimensionalization) . Therefore, the charge background is given by 

(3.1) m e (x) = s~ 3 m (x/e), 

where m a is T-periodic. Note that the factor e~ 3 comes from rescaling so that the 
total background charge in one unit cell is the constant Z independent of e, i.e. 

(3.2) / m £ (x)dx = Z. 

We introduce the following notations for cell average in physical and reciprocal 
spaces 

</(*)>* = I /(*) d*. £ g(k) dk = ±-J^ g(k) dfe. 

When there are no external applied potentials (V ex t = 0, A ext = 0), TDDFT system 
can be written as 

(3.3) i^f- = |(-feV - eA') 2 r 3 + V e ^ e v 

(3.4) - A(j) £ = £(p £ {t,x) -m £ {x)), 

(3.5) ^-A^ + |(V^^ 2 J E , V-A £ =0, 

(3.6) V e (t, x) = (j) e {t, x) + r/( £ y (t, x)). 

We assume that there exists a ground state for the unperturbed system, with 
density having the lattice periodicity 

(3.7) p £ (x) =e- 3 Pgs (x/e), 

where p gs T-periodic. Absence of external perturbation implies that the system will 
stay at the ground state with no electronic current and hence no induced vector 
potential, 

J £ (t,x)=Q, A £ {t,x)=Q. 
The evolution equations are then simplified as 

Bib ■ f 2 

(3.8) i _^ = __ A ^ + ^(x)^. ) 
(3-9) V £ ^cf> £ (x)+r 1 (p gs (x/e)), 
(3.10) -A<j) £ = £ - 2 (p gs (x/e) - m (x/e)). 

Note that the potential is independent of time if there is no external perturba- 
tion. It is easy to see that the potential is er-periodic. We denote the potential 
corresponding to the ground state as V £ (x) — v (x/e) — v gs (x/e) where v gs is 
T-pcriodic. 
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The Hamiltonian operator for the ground state is independent of time, given by 

(3.11) H s =~A + v ss (x/e). 
Define the rescaling operator 5 e as 

(3.12) (SJ)(x) = e-^fix/e). 

It is easy to check that 8 e is a unitary operator. We have 

(3.13) H = -^A + v Q (x) = S* £ H^S £ . 

Since Vq is T-periodic, H is invariant under the translation with respect to the 
lattice L. The standard Bloch-Floquct theory gives the decomposition of H a , 

(3.14) H = I H 0)k dk, 

where i?o,fc is an operator defined on Ll(T) for each feel*, 

L 2 k (T) = {/ £ L 2 (T) | r R f = e~ iR k f, VJ? 6 L}. 

Here tr is the translation operator, i.e. Tnf(x) = f(x + R). The operator i7 ,fc 
has the spectral representation 

(3.15) i? ,fe = ^-En(fc)|V'n,fc}(V'n,fc| ; 

n 

where E n (k) is the n-th eigenvalue of i?o,fcj and tp n ,k is the corresponding eigen- 
function (named as Bloch wave in literature) with 

w„, fe (x) = e- lkx ^ k (x) 

being T-periodic. Moreover, the spectrum spec(-ffo) has the band structure, 

spec(ffo) = IJ U E ^ k )- 
n fcer* 

Denote the spectrum for the first Z bands by <rz, 

z 

(3.16) u z = (J (J E n (k), 

n=i fcer* 

where E n (k) is the n-th eigenvalue of Hq. 

We assume that the ground state satisfies the gap condition, 

(3.17) dist(er z ,spec(iJ ))W) = E g . 

In physical terminology, the system is called a band insulator with band gap E g . 
For convenience, we use the bra and ket notations 

</(C)|£| 5 (OW) = I f*(C)ICg(C) dC, 

where JC : L 2 (T) — > L 2 (T) is a linear operator. 
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3.2. Macroscopic perturbation. We are interested in the dynamics of the elec- 
tronic structure in the presence of the external potentials A ex t(t, x) and V cx t(t, x). 
We assume that A cxt and V cx t are smooth functions in both t and x and periodic in 
space in the domain T. Hence, the characteristic length scales of external applied 
fields are 0(1), while the lattice constant is 0(e). We consider the continuum limit 
e — >• 0; the disparity of the space scales leads to macroscopic Maxwell equations. 

We consider the following system with periodic conditions on T, 

dip? 

(3.18) i-^ = H £ ^, 

(3.19) - A<j) e =e(p £ (t,x) -m (x/e)), 

(3.20) ^A £ — AA E + (V<f ) = e 2 J £ , 

(3.21) V-A £ =0, 

(3.22) V E (t, x) = </) E (t, x) + v(e 3 p E (t, x)), 
where the Hamiltonian operator H e is given by 

H £ = ~(-feV - e{A e + A cxt )) 2 + V £ + V cxt . 

We have used the superscript e to make explicit the dependence on the small 
parameter. The density and current is then given by 

Z/e 3 Z/e 3 
k=l k=l 

Here Z is the number of electrons in one unit cell, which equals to the background 
charge (|3.2p . We remark that in the domain T, since the lattice constant is e, there 
are e~ 3 unit cells in total, and hence TV = Ze~ z electrons under consideration. 

3.3. Main result. Define the limiting macroscopic potentials as 

(3.23) l7o(*,a!) = lim(V £ (t,x) + V cxt (t,x) - v gs (x/e)), 

(3.24) A (t,x) = lim(A E (t 1 x) + A cxt (t,x)); 

e^O 

and the corresponding electric and magnetic fields 

Q 

(3.25) E(t, x) = -V x U (t, x) - -A (t, as), 

at 

(3.26) B(t, x) = V x x Ao(t, x). 
Define the electric field in frequency space, 

y* oo 

E(u,x) = / e iut E(t, x) dt, 
Jo 
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and similarly for B, Uq and Aq. We have 

(3.27) E(u, x) = -V x Uo(oj, x) + iujA (uj, x), 

(3.28) B(w, x) = \7 x x A (uj, x). 

We will show that TDDFT system gives arise to the effective Maxwell system as 

(3.29) V x ■ (S(u)E(u, x)) = p^ t (oj, x), 

(3.30) Va, -B(ui,x) = 0, 

(3.31) V x x E(u),x) = iioB(uj,x), 

(3.32) V x x B(ui, x) = —iuj£(uj)E(Lu, x) + J cxt (w, x), 

where 



Pext(w,a;)= / e iwt pext(t, x)dt, 
Jo 

poo 

J™t(w,aj)= / e ia,t J ext (t,a;)dt ) 
Jo 

with poxt and J oxt given by 

p ext (t,x) = -AajFext^x), 

Jcxt(t,x) = ^A ex t - A^Aoxt + — (VajKxt) • 

The system (|3.29[) - (|3.32p are (nonlocal) Maxwell equations with dynamic dielec- 
tric permittivity matrix £ a/ 3 = 5 a p + A a p given by 

-f — ; 7TT(un.k\idkJu m _ k ) r2(r Ju n _ k \idkg\u m _ k ) L 2, r) dk 

lr* w + u! mn (k) ^ ' p v ' 

n<Zm>Z Ji m " v ; 

^-L Jr* w -aw(fc) ■ w i (i) 

n<Z m > Z 

—3m j- (u nt k\idk a \u mt k)(u ni k\idkJu m ,k)dk 

UJ — — /p* 



Here the potential operator V is the linearized effective potential operator at the 
equilibrium density po : 

(Vf)(z) = cl>(z) + r 1 '(po(z))f(z), 
-A z cf>(z)=f(z), (4>)=0. 
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The operator \u> and the function / are defined as 

Xu>v = ~y2 V f - 1 lu , u nM u* mk {u nM \v\u mM ) L 2 {r) dk 



1 


U! 


+ uj mn (k) 




1 


U) 






1 


UJ 


+ uj mn (k) 




1 







/M = - Y] T 1 ~ 77T u ™.fc M m,fc( u r ! .fc|iVfc|u mi fc) L 2 (r) dfc 

n<Z m>Z 

+ V V f tttU^ k u m . k (u n ^\iV k \u m .k) r2(r] dk. 

T 7 Jr> w ~ w„ m fe ' n v-i 

n<Zm>Z x ' 

Remark that the dynamic permittivity matrix £ is completely determined by the 
linear response of the unperturbed electronic structure. 

The Maxwell equations (|3.29[) - (|3.32[) are nonlocal in time, since the permittivity 
£ depends on uj. While the system is local in space due to the limit e — > 0, the 
nonlocality in time is natural since there is no scale separation in time. 

We also note that while we obtain a nontrivial effective permittivity, the effective 
permeability in the equation equals to 1, the same value as in the vacuum. Phys- 
ically, this is consistent with the situation of semiconductors or insulators under 
consideration here. 

4. Asymptotic analysis of the Schrodinger-Maxwell equations 

To derive the effective Maxwell equation, let us take the following ansatz for the 
system ((3~T8| - ((3~22"]) . 

(4.1) p e (t, x) — e~ 3 p (x/e) + s~ 2 pi(t,x,x/e) + e~ 1 p 2 (t,x,x/e) H , 

(4.2) J £ (t, x) = s~ 3 J Q (x/e) + e~ 2 Ji(t, x, x/e) H , 

(4.3) A e (t, x) + A cxt {t, x) = A (t, x, x/e) + eAi(t, x, x/e) H , 

(4.4) <p £ (t, x) + V cxt (t, x) = <t> (t, x, x/e) + e<j)i(t, x, x/e) 

+ e 2 4>2(t, x, x/e) H , 

where the higher order terms are omitted. We also assume that the dependence on 
the fast variable z — x/e is periodic for all these functions. 

The main strategy of asymptotic analysis is as follows. We first apply a two-scale 
expansion on the Maxwell equations (|3.19[) - (|3.20[) . which produces the asymptotics 
of Hamiltonian; then by Dyson series we obtain the asymptotics of density and 
current; the effective equations in time domain are derived by taking the ^-average 
on the second order perturbation equations of the Coulomb potential and vector 
potential. The asymptotics is somewhat nontrivial. The Coulomb interaction makes 
the leading order potential dependent on the macroscopic average of the third order 
density. To close the asymptotics, one has to show that the macroscopic average 
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of the third order density only depends on the leading order potential, but not 
on higher order terms of the potential. This amounts to establishing the local 
neutrality of the system, which will be explained in detail below. Finally Fourier 
transform gives the effective Maxwell equation in frequency domain. Note that we 
have assumed that the leading order density and current only depend on the fast 
variable x/e. This will be justified by the asymptotics. 

4.1. Asymptotics of the Hamiltonian. For the Coulomb potential, substituting 
the ansatz (|4.1j) and (|4.4|) in (|3.19j) and organizing the results in orders, one has 

(4.5) - A z (f) Q = p - m , 

(4.6) - A z 0x -2V X • V Z O = pi, 

(4.7) - A z <f> 2 - 2V X • V z 0i - A x <f> = p 2 + p cxt - 

Recall that p cxt {t,x) = -A x V cxt (t,x). 

For the exchange-correlation potential, Taylor expansion yields 

r](e 3 p e ) = r){p ) + er)'(p ) Pl + \s 2 r)"(p )pl + e 2 i 1 '{p )p 2 + ■■■ 

( 4 - 8 ) 

= vo + em + £ I ^2^ — , 

where the last equality gives the definition of rji(t,x,z), 

T) (z) = 77(A)(2)), 

771 (i, x, z) = 7}'(p (z))pi(t, x, z), 

r) 2 (t,x,z) = ^i 1 "(p (z))p 1 (t,x,z) 2 +r)'(p a (z))p 2 (t,x,z), 

and similarly for higher order terms, which we omitted in the expression. 
Therefore, the total potential V s can be written as 

V e = c/> + V ext + V 

(4.9) = (<£o + Vo) + e(0i + m) + £ 2 (^2 +»») + ■■• 
= V Q + eV x + e 2 V 2 + --- . 

Similarly, we write down the equations for the vector potential using ansatz 
and JO}: 

(4.10) - A Z A = 0, 

(4.11) - A Z A X - 2V X ■ V Z A + ^(V z Vo) = J , 
d 2 

(4.12) — A - A X A Q - V Z A 2 - 2V X • ^ Z A X 

d 

+ ^(V z Vi + V.Vfc) = J 1 + Joxt- 
Recall that J ext (t,x) = ^zA cxt - A x A cxt + ^ (V x V C xt)- 
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Note that the solvability conditions of (|4.5p - (|4.6|) and (|4.11j) impose the following 
constraint on po, pi and Jo, 

(4.13) (p {z)) z = {m (z)) z , ( Pl (t,x,z)) z = 0, (J (*)>» = 0. 

Therefore, the Hamiltonian operator can be written as 

H £ (t) = -±e 2 A x + v (x/e) + U (t, x) + evi(t, x,x/e) 
u 14 n +eU 1 (t,x) + ie 2 A Q (t,x) ■ V x + e 2 V 2 (t, x, x/e) 

+ e 2 lMt 2 X)l2 + ze'A^ x, x/e) ■ V x + • • • 

where we omit the higher order terms. In (I4.14[) . we define vq, v\ and Uq, U± as 
the microscopic and macroscopic components of Vq, V\ respectively, 

U {t, x) = (V (t, x, -)) z , v (t, x, z) = V (t, x, z) - U (t, x); 

(4.15) 

U x (t, x) = (V 1 (t,x,-)) z , Vl (t, x, z) = Vx (t, x, z) - U x (t, x). 

4.2. Asymptotics of the density and current. The initial state is given by the 
ground state of the unperturbed system, which implies that the density matrix is 
given by the projection operator to the occupied spectrum of the ground states in 
the beginning (with lattice parameter e), 

(4.16) p e {0) = V e . 
Then, the density matrix at time t is given bj@ 

(4.17) p e {t) = Tcxp^-i^ H E {T)dT^V £ (^rcxp(^-i £T(r)dT^ , 

where T is the time ordering operator. Therefore the density is given by the 
diagonal of the kernel of the operator p e , 

(4.18) p E (t,x) = p £ (t,x,x) = Tcxp(^-iJ ir(r)dT^ 

xP E (Texp(^-i H E (r)dT^ (x,x), 

where the right hand side means the diagonal of the kernel associated with the 
operator. The ground state density is given by p^x) = P E (x,x). 

We investigate the asymptotic expansion of p E (t,x) determined by the Hamil- 
tonian (|4.14[) . The equation (|4.18[) implies, to obtain the density at (t,x), we first 
evolve the system backwards to the initial time, project it onto the ground states 
of the unperturbed system, then evolve it forwards in time to t at x. Since the 
current time scaling is 0(1), when e goes to 0, the domain of dependence and do- 
main of influence are also of the scale 0(e) for the system in the time evolution. 
In other words, the density at (t, x) only depends on the Hamiltonian of a small 



2 In the language of physics, we are using the Heisenberg picture. 
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neighborhood (0,i) x B(x,0(e)), where B(x,0(e)) indicates a ball centered at x 
with the radius O(e). 

Accordingly for a fixed point x e T, we expand the Hamiltonian operator H 
around x. For clarity, we write H as an operator on L y (T), 

H e (t) = -±e 2 A y + v (y/e) + £/ (i, y) + ewi (i, y, y/e) 

+ eU x {t, y) + ie 2 A {t, y) • V y + e 2 V 2 (t, y, y/e) 

+ e 2 lMt 2 y)l2 +ie*A 1 (t, y ,y/e).V y + --- . 

Let i?o(*> 33 ) be the leading order operator 

(4.19) floras) = -±e 2 A y + v (y/e) + U Q {t,x). 

Here H^t.x) is an operator on L 2 with (i, x) as parameters. Similar notations 
will be used throughout the paper. Denote the difference as 5H £ (t,x) = H £ (t) — 

5H £ {t, x) = U (t, y) - U (t, x) + evi(t, y, y/e) + eUi(t, y) 

( 4 - 20 ) + ie 2 A {t, y)-V v + e 2 V 2 (t, y, y/e) + e 2 ^ 

+ ie 3 A 1 (t,y,y/e)-\7 y + -- - . 

Expand 5H e (t,x) into orders, and the first two orders are 

(4.21) SH e {t,x) = {y-x)-V x U {t,x)+ev 1 (t,x,y/e)+eU 1 (t,x)+ieA (t,x)-e\7 y , 



and 



5H £ {t,x) =\{{y - x) ■ V x ) 2 U Q (t,x) +e(y - x) ■ V x V!(t,x,y/e) 
( 4 - 22 ) + e(y x) ■ VMt, x) + e 2 V 2 (t, x, y/e) + £ 2 ^f^ 



i(e(y - x) ■ V x A (t,x) + e 2 Ai(t,x,y/e)) ■ eV y . 



Recall that we are approximating the Hamiltonian in a ball around x with the 
radius O(e), and hence (y — x) is treated as 0(e) in the above expansion. 



14 
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We expand the time propagation operator as 



Texp^-z^ H £ (T)dT^j =U^(x) 



-ij^ Texp^-i J H £ (s)ds^SH £ (T,x)U^(x)dT 

=U*$(x)-i f t U^(x)SH £ (T,x)Uf Q (x)dr 
Jo 

- f r 2 Uli(x)5H £ (r 2 ,x)U £ ; Tl (x) 
Jo Jo 



'0 Jo 

x8H e {T U x)U s T f Q {x)dTid.T 2 



+ 



where U^° t2 (x) is the propagation operator corresponding to H e , 

U e t f M (x) = T exp (-i jf 1 H%(t, x) dr j . 

Therefore by (l42l"l) - (|4~22)) . 

Texp^-ij^ H E (r)\ dT=U^(x)-i U £ t '°{x)5Hl{T,x)U e T fl {x)dr 

- i f Uf T (x)SHI (r, x)Uf (x) dr 
Jo 

(4.23) r t r r 2 

/ / U^ 2 (x)SHt(r 2 ,x)U^ Tl (x) 
Jo Jo 

xSH £ 1 (T 1 ,x)U E T '° (x)dT 1 dT 2 + --- 
=U$(x) + eU^ix) + e 2 U^{x) + ■■-, 

where the last equality defines UIq(x) for j — 0,1,2. Higher order terms can be 
written down similarly. 

Substituting the expansion (|4.23|) in (|4.18[) . we obtain to the leading order 

p s (t,x) = (u$(x)V e U$(x))(x,x), 

and also to the higher orders, 

P \(t,x) = (hl^(x)V e U^(x))(x,x) + (u^(x)V s U^(x)y X) x), 

p%{t,x) = (hl^(x)V e U^(x))(x,x) + (u$(x)V s Ul : l(x))(x,x) 

+ (l($(x)V e U$(x))(x,x). 

Therefore the density is given by 

p e (t, x) = p s (t, x) + split, x) + e 2 p\{t 1 x) + 0{e 3 ). 
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The current is given by 



1 



J e (t,x) = -[eV y , p s (t)] + (x,x), 

where p £ (t) is the density matrix in (14. 17)) and [a, b] + = ab + ba is the anticomrnu- 
tator. 

With the help of (14.231) . similarly one has 

J £ (t,x) = J%{t,x) + eJ\{t,x) + 0{e 2 ), 



where 



J °=2* 
1 



eV y , U s t ${x)V*U${x) 



(x,x), 



2i 



eV y , U^{x)V e Ul : t{x) + U e t $ {x)V e U^{x) 



6,1 1 



(x,x) - plA (t,x). 



4.3. Rescalings of Hamiltonian, density and current. Notice that the oper- 
ator H^{t, x) in (|4.19p can be rescaled as 

8 e H$(t,x)S* e =H (t,x), 

with 5 e as the dilation operator and Hq is given by 

H (t,x) = -|A C + «o(C) + V (t, x), 

where C = vl e is the small scale spatial variable. 

Therefore, we can rescale the expressions of p £ (t, x) and J^(t, x) by 

Pj(t,x) = e~ 3 pj(t,x,x/e), Jj(t,x) = e^ 3 J 1 {t,x,x/e), j = 0, 1, • • • , 

and 

(4.24) Po (t, x, z) = (l$ >0 { X )VUl t {x)) (z, z), 

(4.25) Pl (t,x,z) = (ul Q (x)VUl t (x,z)){z,z) + (ul (x,z)VU^ t (x))(z,z), 

(4.26) P2 (t, x, z) = (uUx)VUl t (x, z)) (z, z) + (ul (x, z)VU^ t (x, z)) (z, z) 

+ (ul (x,z)VUg >t (x))(z,z), 

V c , UUx)VUL(x) 



(4.27) J (t,x,z) 

(4.28) Ji(t,x,z) 



1 

~2~i 
1 



(z,z), 



(z,z) 



V C , Ul G (x)TUl t {x,z)+U} fl {x,z)VUl t {x) 
- p (t,x,z)A Q (t,x), 

where V is the rescaled (lattice parameter 1) density matrix for the unperturbed 

2 

c 

U° i3 (x) = T 'exp(-i / ffo(r,a;)dr), 

J S 



system, and IA\ s are propagation operators defined on L\ by 
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(4.29) Ul s (x,z) = -if UI t { X )5H x {t, X ,z)U%{x)At, 

J S 

(4.30) 

Ul s (x,z) =-iJ U° :T (x)SH 2 (t,x,z)U%(x) dr 

- j j 2 Ul T2 (x)SH 1 (T2,x,z)U° 2tTl (x)SH 1 (T 1 ,x,z)U° uQ (x)dr 1 dT2, 

in which 

(4.31) H (t, x) = -±A C + v (C) + U (t, x), 

(4.32) 5Hx(t, x, z) = {C-z)- \7 x U (t, x) + vi(t, x, C) + Ui{t, x) + iA ■ V c , 

(4.33) 5H 2 (t, x, z) = \ ((C - *) • V„) V (t, x) 

+ (C - z) ■ V,(wi(t, x, C) + Utit, x)) + V 2 (t, x, C) 

+ i((£-z) ■ V X A Q + A^j -V c + i|A | 2 . 

Notice that the dependence of the operator Hq(t,x) on x lies only in U(t,x) 
which works like a number as an operator on L 2 ^. It implies 

Ut s (x) — exp ^— i J U(t, x)^j exp(—i(t — s)H ), 
where Hq agrees with the unperturbed Hamiltonian 

(4.34) H = -±A c +v (O- 

Moreover, in the expression of pj, the phase factor exp(— i J U(t, x)) will not ap- 
pear since it gets canceled by its complex conjugate. Therefore, we can simply 
take 

(4.35) W^expHt-s)^), 

which is, in particular, independent of x. 

It results that the leading order density agrees with the ground state electron 
density of the unperturbed system, 

(4.36) Po(t,x,z) = e- UH «Ve ltHo (z,z) = p gs (z). 
This also shows the following proposition. 
Proposition 4.1. po is independent oft and x and 

(po)z = (m ) z , 

which satisfies the first constraint of (|4.13[) self- consistently. 

5. Effective equations in the time domain 

In this section, we derive the main result in time domain. It will be connected 
to the main result described in frequency domain in the next Section. 
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5.1. Effective equations in time domain. Let us first summarize the resulting 
equations in the time domain. Recall that 

V (t,x,z) = v (z) + U (t,x), 
Vi(t,x,z) = vi(t,x,z) + Ui(t,x), 
A (t,x,z) = A (t,x), 

then one has the following effective equations from (|3.18p - (|3.22[) . 

• The microscopic scalar potential Vq(z) is the same as that of the unper- 
turbed system. 

• The potentials v%(t,x,z), Uo(t,x) and Ao(t,x) form a closed system as 
described below. The microscopic scalar potential i>i(t, x,z) is given by 

ui = Vi(t,x,z) - (Vi(t,x,z)) z , 

(5.1) Vi = Vpi = 4>i +r)'{po)pi, 

- A z 0i = pi. 

Recall that V is the linearization of the effective potential operator at equi- 
librium density. Here, the density pi is given by the equation 

(5.2) (I- X V) P1 = [ f(t-T)-V se U (T)dT+ [ g(t-r) ■ A (r)dr, 

Jo Jo 

where / is the identity operator, 
XVpi = X vi = -23m E E / e^" (fc)(t - T) 

n<Z m>Z'' •' T * 

x u n ^u* n k (u n . k \vi(T)\u mik ) L 2^) dfcdr, 

and 

f{s) = -23m E E f ei " m " (fe)SM n,fc<n,/=( u ™,feN V fcl u m,fe)L2(r)dfc, 

n<Z m>Z 

g(s) = -23m E E T ei " m " (fc)Su ™,fc M m,fc( u »,fel* v cl u ™,fe)£ 2 (r) dfe. 

n<Z m>Z •' F " 

Here we have introduced the short hand notation 
Wmn(fe) = E m (k) - E n (k). 

• The macroscopic scalar potential Uo(t,x) satisfies 

(5.3) - AU - P a pd Xa d Xf> U = Q a {d Xa vi) + R a p(d Xa (A )i3) + p ext (t,x). 
Here 

P*p (d^d^Uo) = 23m E E / / e <Wm " (fc)(t - T) (u„, fc |ia fc J« m , fc > i2(r) 

X (Un,fc|iQfe J a|lim,*}ia(r)ft«>a^B)3^Q( r ) dfc dr ' 
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a (d Xa v 1 )=23mY / Y, I f eiUmn{m ~ T) (un,k\idkJurn,k) LHr) 

n<Zm>Z^° 

x {u n ^\d Xa vi(T)\u m ,k)L 2 (T) dfcdr, 



and 



R a0 (d Xa (A o ) )=23mY, J2 I { z i " mn(k){t ~ T) K,k\idk a \u m , k ) 1 

n<Z m>Z J° 



ILHT) 

x (w„,fc|i9 C/3 \u m . k )L^{T)d Xa {A )p dfc dr. 
• The macroscopic vector potential A$ satisfies 

(5.4) — (A ) a -A x (A Q ) a + —(d Xa U ) = 

S a (vi) + M a p(d Xl3 U Q ) + N a p((A )p) - (po)z(A ) a + {J cyLt ) a {t,x), 
V x • A Q = 0. 

Here 

S a ( Vl ) = 23mJ2J2 l{ eiWmn(fe)( * _T) K^RK^) i2( r) 

n<Zm>Z'' 

x (u7i,fc|«x(T)|«m,fc)i 2 (r) dfcdr, 

M^^Db) =23fm ^ J2 f { e iUm " (fc)(t " T) <"„, fc |i3cJ u m,fc> La(r) 

x (M„,fe|«9 fe/3 |um,fc)L2(r)5x^t / "o('r) dfc dr, 



and 

JV a /»((A>)/») = 23m J] E / / e --»W(*-)^ nife |0 c J Umife ) i3(| | 



n<Z m>Z 



x (u„. k \id^ |u ro ,fe)z,2(r)(>*-o)^(r) dfc dr. 

The above equations from <|5. 1|) to (|5.4[) form a closed system that determines the 
macroscopic potentials Uo(t,x) and A (t, x). 



5.2. Derivation. We first give a description on how to derive the equations (|5.f P ~ 
()5.4|) . By taking the ^-average of (|4.7[) and ()4.12|) . one can get the effective equations 
for Uo{t,x) and Ao(t,x). However, in order to get equations in explicit form, one 
needs the expressions of pi, {pi)zi Jo an d (Ji)z in terms of Bloch wave function 
{ipn,k} or its periodic part {it nj fc}. This will require the following three lemmas. 
The first lemma states the property of perturbed density under Hamiltonian per- 
turbation, and the other two introduce some identities related with Bloch waves. 
These identities will be useful in simplifying the expressions of density and current. 
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Lemma 5.1. We consider the electron dynamics under the perturbed Hamiltonian 
H = Ho + Y]j—i ''{Vj + iAj—i ■ V$) and denote the density as 

p(t,x,z) = p (t,x,z) + ^2e k p k (t,x,z). 



Assume the initial condition is 

p(0,x,z) = p ss (z), 

then one has for any k, 

(p k (t,x,z)) z = 0. 

Proof. Since we choose the Coulomb gauge so that • Aj-i = 0, j = 1, • • ■ , J, 
the operator 7~exp(— i f Q H(t) dr) is unitary, which produces 

(p(t,x,z)) s = (p(0,x,z)) z = (p gs {z)) z . 

Therefore by (|4.36|) and 

(p(t,x,z)) z = (p (t,x,z)) z + ^2e k {p k {t,x,z)) z , 

k 

one gets, for any k, 

(p k (t,x,z)) z = 0. 

□ 

Lemma 5.2. For any n, m 6 Z + and k,p £ T* , the following equations hold in 
the distributional sense. 

(5.5) / (C - ^OC/cVVp d C = \ 5 nm {~z - id p ) 

+ (un,fc|iV fc |u mi fe} L 2 (r )|(5(p - fe) |r*| , 

(5.6) / vi(r, a;,C)C fcVWpdC = {u n , k \vi(T)\u m . k ) L 2 {r) S(p - fe) |r*| , 

( 5 - 7 ) / i>n fc(*V«;)'!/'m,pdC = (u„,fc|zV C |u m ,fc) L 2 (r) (5(p - fe) |r*| . 

JW.3 

Proof. By the definition of Bloch wave, direct calculation yields 
/ (C-*Xh^ n HpdC= / (C-^)e< ( f- fe )<. fcUm , p dC 

= £ /(C + X.-^+^H^^u^dC 



x.<eh 



20 WEINAN E, JIANFENG LU, AND XU YANG 

For the right hand side, we have 



(5.8) 



[ e^+ x ^-^u* ntk (id p )u m!P dC 
<r /-tt ** r 



X,-eL 



= -^pf E e^ (p - fe) f e< (P- fc )< fc u m , p dc' 



To further simplify the above expression, we use the Poisson summation formula 
(5.9) E e lX ^ = |r*| E *( fe -^)> 



x^gl 



in the distributional sense. Substitute (|5.9j) into f|5 . 8[) . we obtain 
E / (-id P e i ^ x ^- k Au* n , hUm , p dC = 

ir. c irVr\ / 

- |T*| zd p ^(p - fe) jf e^ ^- fe )< ife u m , p dc) 
+ |r*| 5(p -k)J e< ^u* nM {id p )u m , p &C 

- |r*| id p Up -k)J e< (p- fe )< :fcUm , p dC 

+ |r*| S(p - fe)(M„. fe |iV fe |u TO , fc ) L 2 (r) 
= - |r*| id p ( S(p - fe)(itn,fc|«m,fe)r,2(r) 



+ |r*| <5(p - fc)(u„, fc |zV fe |u m , fe ) L 2 (r ) 
= |r*| (-id p )S(p - k)5 mn + \T*\S(p - k)(u n! k\iVk\u m ,k)^(r)- 

The last equality follows from the orthogonality of {u Ui k} for each k. 
Similarly we have 

E / -^w^wdc 

(5.10) = -a |r*| 5(p - k) J e«P-Vu* n k u m . p &t 

= -Z |r*| 8{p - k)(u n> k\u m ,k)L2(r) 

= -z\T*\S(p-k)S mn , 



EFFECTIVE MAXWELL EQUATIONS FROM TDDFT 



21 



and 

/ v 1 (T,x,C)ipZ k ip m!P d£= / i;i(r,a;,C)e iC ' (p ~ fc) ?v fc u m ,pd< 

Jr 3 Jk 3 

= (u„,k\vi(T)\u m , k ) L 2( r )6(p - fc) |T*| , 

where we have used the periodicity of i>i(t, x, £) in £. 

Hence combining the above equations together yields (|5.5[) , and the last equality 
proves (|5.6|) . The proof of (|5.7|) is essentially the same as (|5.6|) which we will omit 
here. □ 

Lemma 5.3. For any n,m G Z+ cmc? k,p G T* 7 we /iawe in the distributional 
sense, 

/ (C - Z)«(C - Z)p1pn,k1>m,p dC = 
JR 3 

*nm(Zo2/i - d Pa d Pli )S(p - fe) |T*| 

( 5 ' 11 ) + (■"n,fe|i9fe (3 |it m> fe) z ,2( r) (-i2; Q , - id p J5(p - fe) |T*| 

+ (w„,feK9fe„|M m ,fe)i=(r)(-*2;/9 ~ id Pfl )6(p - fe) |r*| 
- (wn,fc|9fc Q 9fc3|Mm,fc)L2(r)<5(p-fc) |r*| . 



/ 



(C - z)^^ ui(t, a;, OVCfeVw dC = 
( 5-12 ) (-" n , fc |5 XQ t;ii5fc Q |u m ^) L 2( r) (5(p-fc) |T* 

+ {u n .k\d Xa v 1 \u m . k ) L 2 {T) (-z a - id p JS(p - fc) |r* 



(5 - 13) K >fe |i9 fet ,(-fe^ + i9 C/3 )|u m>fe ) L2(r) 5(p-fe)|r*| 

+ (un.klid^^Um^) L 2 (r) (-z a - id Pa )S(p - fc) |r* 

Proof. We calculate 

/ CaCsC.fcVW dC = / CC/3< l ,feMm,pe i(P " fc)C dC 
JR 3 JR 3 

= - E I <^ P d P j P ^- k ^+ x ^ dc 

V /—IT T 
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The Leibniz rule gives 
Jr 3 

s ^(/< fe w i(p - fe)c dcE ei(p " fe)x 



9 Pa 



r -V- *— t 



Therefore applying the Poisson summation formula (|5 .9[) yields 

/ CaC^CfeVW d C = - S nm d p d p S(p-k) \T*\ 

+ (Un,k,d kf! U m ,k) Li(r)dpJ{p - k) \T*\ 
+ (un,k,dk a u mi k) L 2( r )d Pfj d(p - fc) |r* 
- (u n . k ,d ka d kfl u mik ) L 2 (r) 6(p - k) |r*| , 

and hence using (I5.5[) and (I5.10[) . we have (15. lip . 

The calculations for (|5.12l) and f|5.13[) are analogous and omitted here. 



□ 



5.2.1. Derivation of the equation ([Oj) , By (|425|) . (|435j) and (|4T29|) the first order 
density perturbation reads as 

ft - (^ i0 (x)7^x, *))(*,*) + (ul (x,z)VUS, t (x))(z,z) 

(5.14) t 

fop / e* rHo *Fi(T,x,*)e i ( t-T ) fl ' dT(*,*)+c.c., 

J 



where we have used the fact that 



(l4 t0 (x,z)-PU$ it (x))(z,z) = (U° fi (x)VUl t (x,z))(z,z), 
as a direct consequence of 

(Ul (x, z)VUl t {x)Y = U° (x)VUl t (x, z) 
in the operator sense. 
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Substitute (|4.32j) into (|5.14j) . we obtain 
(5.15) 

pi(t,x,z) =-23me- itHo V [ e lTli °{C - z) ■ V x [/ (t, a;)e l(t ~ r)HQ dr(z,z) 

Jo 

-23me- UHn V f e iTff °«i(r, x, C,) e l{t ~ T)H ° dr(z, z), 
Jo 

-2Jme- HHo V f e lTHo A a (T lX ) ■ ({V c y (t -^ Ho dr{z,z), 
Jo 

in getting which, we have used 

-23me- itHo V f e iTB °Ux{r, x)e t(t - T)H ° dr(z, z) 
Jo 

= -23me- ltHo Ve UHo (z,z) f U 1 (T,x)dr 

Jo 

= - 23mP [ U 1 (T,x)dr = 0. 
Jo 

The first equality above follows from the fact that Ui(t, x) is a number as an 
operator on L^. 

Proposition 5.4. The average of p\ with respect to the microscopic scale vanishes, 

(5.16) ( Pl (t,x,z)) z =0, 

which satisfies the second constraint of (|4.13[) self- consistently. 

Proof. In (|5.15[) the first term of the right hand side is an odd function in z, hence 
when taken the average over z, it gives zero. The second term is the first order 
density perturbation p\ if one takes J = 1 in Lemma 15. 1[ whose average over z is 
also zero. 

□ 

For a more explicit expression for p%, we substitute the spectral representation 
of operator H into (|5.15p . 

(5.17) H o=Y,-f 
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This gives 
(5.18) 

pi(t,x,z) = 



23m E E f / ^wc,p(^ i(E ™- E "' fe,(H 

x f (C-zW nk (C)^ m AC) d Cdkdp-V x U Q (T,x)dT 

JjgL 3 



23m EE/*/ ^, fe (^)c P (^)e i(E "" p - s - fe)(i - T) 



x / w 1 (r 1 x 1 C)Cfc(C)V'm,p(C)dCdA!dpdr 3 

^ ee/7 ^(*,p(^"" mm 

' — /o J(r*) 2 



n<Z m>Z 



/ < fc (C)(iV c )V m ,p(C) • Ao^^dCdfcdpdr. 

JR 3 



Substituting §5J^-$5J} in ([5T8|) . we obtain 



t 

* . I -\ p iu mn (k)(t-T) 



n<Z m>Z 



JT* 

x (w n! fc|iVfc|w mifc ) L 2 (r ) dfc • V x U (t, x)dr 
t 

* . I y\p l ^ m -n(k)(t-T) 



- 23m V V / f u n , k {z)u* mtk {z)e 
(5.19) n<z m >z Jo Jr ' 



x (u nt k\vi{T, xX)\u m ,k) lz(t) dfcdr, 
23m EE/*/ u n , fe (z)< fe (z)e^»«('-^ 



x (w nifc |iV c |u m! fc) L 2 (r ) • A (t,x) dkdr. 



This implies (f572|) . 

Remark. From (j5.19[) and using the orthogonality of {u n ^k} for each fc, we once 
again see that 

(pi(t,x,z)) z = 0. 
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5.2.2. Derivation of the equation ([Oil . By (j4T26|) . (|435|) and (jOO]) . the second 
order density perturbation reads as 



p 2 {t,x,z) =-2V\ze- UHo V [ f 2 e lTlHo SH 1 (T ll x,z)e^ T2 - T ^ Ho 

Jo Jo 

x SHxfax, z)e i{t - T * )Ho dn dr 2 (z, z) 

-23me- i * ffo 7' / e lTHo 6H 2 (T, x, z ) e l ^ T ) Ho dr(z, z) 
./o 

+ / e- l (*- Tl ) H °,5 J ff 1 (r 1 ,a;,2)e- iTlHo dr 1 7' 
Jo 

Jo 



X 



Since commutates with exp(-itHo), we may simplify the above expression as 



p 2 (t,x,z) = 

-29\eV f f 2 e^ Tl - t)Ho 6H 1 {T 1 ,x,z)e l(T2 - Tl)Ho 
Jo Jo 



(5.20) 



x SH^x, Z ) e ^-^ H " dn dr 2 (2, z) 

-23m7> f e l{T - t)Ha SH 2 (T,x,z)e lit - T)Hn dT(z,z) 
Jo 

+ f e- l{t - T ^ Ho SH 1 (T 1 ,x,z)P 
Jo 

x f e l{T2 - Tl)Ho 5H 1 {T 2 ,x,z)e l ^- T ^ Hn dr 2 dri(z,z). 
Jo 



Proposition 5.5. The average of p 2 (t,x, z) is given by 
(5.21) 
{p 2 {t,x,z)) z = 

rt 



23m (r J e^ H ^((C - z) ■V x ) 2 Uo(T,x)e^ t - T ^ d T (z,z) 
23m (v e^ T ~ t)Hn (C - z) - VvV^xXV^" dr(z,z)^ , 
23m (v J e ^-^ H ° f « - z) ■ V x A (t, x) \ ■ (iV c y ^ H ° dr(z, z) 
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Proof. Substituting the expressions of 5 Hi and SH 2 (|4.32|) - (|4.33|) into (|5.20|) and 
taking average with respect to z, one has 

(p 2 (t,x,z)) z = h+h 

- 23m (v J* e t{T -^ Ho \ ((C - z) ■ V x ) 2 U {t, x)e^- T ^ dr(z, z^j 

- 23m (v e i(T-*)Ro(£ _ z ) . S7 xVi (t, x, Q e <t-r)Ba dr ( Z) 

- 23m (t> j e ^ T -^ Ho ((C - z) ■ V x A {t, x)j ■ (iV ' C )e l(t - T)H ° dr(z, z) 
where I± and I 2 are given by 



x (C-z) • \7 x U (T 2 y H ^- T ^ d Tl dT 2 (z,z) ^ 
x jf e^ota-n)^ _ z ) . V.Ubfa)^' - ^ dr 2 d n (z, z)^ , 



J 2 = -2fHe ( P / / 2 e lHa{Tl - t] {C - z) ■ V x C/o(r 1 )e jHo(T2 - ri) 
Jo 

x L(t 2 ) + A (t 2 ) • (iV^e < ^( i -^dTidr 2 (z,z) 

t /-T2 



x (C - z) ■ V.E/bfo)^*-^ d n dr 2 (z, z)J 
+2^e / jf e **° - 2 ) • V x f/ (Ti)P 



Remark that when calculating (p 2 (t, x, z)) z , we have dropped out the z-average 
of the odd functions in z and the second order density perturbation functions of 
H = H Q + e(v x + Ui+ iA ■ V c ) + e 2 {V 2 + \ |A | 2 + %A X ■ V c ) by Lemma0 
We complete the proof by showing that I x — and I 2 = 0. 
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We denote the second term in I\ as I\$ = $Ke (J* J* K(ti,T2) dn dT 2 (z, z)) z , 
then 

Ji, 2 =9te^ jf 1 X(r 1 ,T 2 )dT 2 dT 1 ( 2; ^)^ 

+ $He^ £ 2 K(T 1 ,T 2 )dr 1 dT 2 (z,z)^ 

= 29U^j[ ^ 2 A"(ri,T 2 )dTidr 2 (z,z)^ , 

where the last equality is obtained by switching n ■<-> r 2 in the first term of 2i i2 
and using the fact that if(n,T 2 ) = if(r 2 ,Ti). Therefore n could be rewritten as 

x (C - z)^ ^^ dn dr 2 (z, z) 

(5.22) t 

+29*e \J J' e iHo(n -^(C - z) a Ve iHo ^ 2 - Tl) 

x (C - z^e^-^L^ dn dr 2 (z, z)\ , 

where we have used the short hand notation U^f T2 = d Xa lIo(Ti)d Xl3 Uo(T2) ■ 
Substituting spectral representation of Ho (15 . 1 7[) into (I5.22j) gives 

h = -2<Ke E ( f r f M'WWe^^K.Hrna, 

n<Z rnl \ J ° J ° "V) 3 

x e^'^^F^^e^^^U^, dn dr 2 dfedpdq\ 

/ Z 

+2<Kc £ E ( f r / ^.p(- z )#, 9 ( 2; ) e ' Bm - p(Tl_t)i? m, P; „, fe 

x e lB "- fc(r2 - Tl) ^ fe ., qe ' iB ^ ( *- r2) C/ T f r2 dn dr 2 dfc dpdq^ , 

where F« k . mp = J Ra (C ~ z) a ^ k (C)^m, P (C) dC, and F" lp . nk , F^ p . tq , F^ k . tq are 
defined similarly. 
We denote I\ as 

h = -21He y, E ^ + 2% EE N X> 

n<Z mi n<Z ml 

then it is easy to see that, by switching (m,p) and (n, k), 

E E E^' +*« E E E^ = °- 

n<Zm<Z £ n<Zm<Z I 
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Therefore 



h = -2*e EEE N %£ + m * E E E N X- 

n<Zm>Z ( n<Zm>Z I 



Making use of the identity (|5.5p produces 

F n,k- m ,p = ^nm(-Z a - id p J + (u», k \idk a \u m ,k) L* (T) j KP - k ) \?*\ , 

f m, p; £, g = |<W(-*;8 -idq ff ) + (u m , P \id Pfj \ut p ) L 2 (r) ^5(q - p) \T*\ . 

Then by the orthogonality of {u Ut k} for each k and using integration by parts 
for the variable p, one could rewrite I\ = + /{ which are given by 



4 1] 



x e ^ m .-(- 2 — O^C-^^ - id q0 )8{q - ky E ^-^U^ T2 dn dr 2 dfcdq) 
+ He EEE(/ / '/ V' m ,fe^%e^- fc ( Tl -*)( Mm , fc | l 9 fc J U „. fc ) L 2 (r) 



x e^to-^M-S/J - i^)<% - k)e lE ^-^U^ T2 dn dr 2 dfcdq^ 
+ 29=te^ E ( / / / ( M ™,feN 5 fc Q l u ™,fe)i 2 (r)(wr i .fc|«9fc^|u m ,(e) 

n<Zm>Z WO ^ Jr * 



L=(r) 

x e i(B mit -B„. fc )(t-n) U ^ 9 dn dT2 dfe \ 

/ 2: 

2^He E E \ / / / ( w m,fcN 5 fc Q l u «.fe)L 2 (r)("m,fe|i9 fe|3 |w n! fc) i2fr > 
^r^V Wo Jo Jr* 



e HBn,u-E m , h )(.t-r 1 ) U ap d Ti dr 2 dfc ) , 
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/{ 2) = -25He V V W / / 2 / Vn,feV , ! fe e i ' B ''.'=( T1 - t )( Un , fe K9 fe Jn ro , fe ) i 2 Cr) 



x e 



iEm,fe(T2— Tl 



+2^e J2J2J2(f I 2 / V'm.fcV'; fc e l£; '"- fc ( Tl - t )(u m , fe |0 fc Jw„, fc ) L 2 (r) 
x e 4f! - fc(T2 - Tl) K, fc N9 fc J^ fe ) L2(r) e^^(*- T2 ){/ T f T2 dn dr 2 dfc 

ft r r 2 



= -2<Ke V] V] \ / / / ("™.fcN9 fcci |u mife ) i 2( r) (u m , fc |z9 fc)3 |u„ ifc ) 

„<Z m >Z W° ^ Jr * 
x e i(B m , fe -E„, i! )(r 2 -T 1 ) Lr a^ ^ ^ dfc ^ 

+2<Ke X! X! ( / / / ( M ™-fel^fc a l u ".fc)i 2 (r)( u r 1 ,fc|«9 fc)3 |u m .fc) L2 , r) 
n<Zm>Z w° J o 

g.tE^-fi^O^-n)^ dn dfc 
By observing that 



we get J-^ = and l[ 2 ^ — since one has the same real part as its complex 
conjugate. 

Therefore I\ = 0. Similar arguments will show that I2 = by making use of the 
identity (|5.6|) . and we omit its details here. □ 



Substituting the spectral representation of H (|5.1T[) into (|5.21|) gives 
(p 2 (t,x,z)) z = 



X 2 



I [ (C-*)o(C-*)^Cfe(C)V'm,pdCdfcdp9 <Ba a a ,pl7o(r,a;)dT 
* it 3 



ri<Z m>Z 



/ (C - *)a0 a!o Ui(' r ' a5 :C)VC JOVVpdCdfcdpdr ) 

JK3 / , 
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By (|5.1ip and (|5.12p and using the integration by parts with respect to p, we 
could simplify the above equality to be 

(p 2 (t,x,z)) z = 

23m £ E f f e^WW 

n<Zm>Z^° •' r " 



x (u n ,k\idk a \ \idk p \u m ,k) L^(T)d Xa d Xl3 U (T) dfc dr 

n<Zm>Z JQ JT " 



ft 



x (u n ,k\idk a \u m ^) L2{r) {u n ^ k \d Xa vi(T)\u m ,k) L^(r) dfcdr, 



+ 2 ^EE f 1 



e tw m „(fe)(t--r) 



n<Z m>Z 



x (M„,fc|i3fc Q |u m: fc) i 2( r) (wn : fc|i3 <)3 |u I „ :fc ) L 2( r )9 a , ci (Ao)^ dfcdr. 
Therefore taking the ^-average of (|4.7[) produces 



- S al3 d Xa d Xll U a = (p 2 (t,x,z)) z + p cxt (t,x) 

= P a p{d Xa d Xf} U Q ) + Q a (d Xa vi) + R a p{d Xa (A )p + Pext(*,aO- 

This proves (J5T 



5.2.3. Derivation of the equation (|5.4p . Substituting the spectral representation of 
H ([5~T7| into (|Q7|| yields 

(5.23) J o=E / ^^n^V^^n.fcdfc. 

n<z ^r* 

Proposition 5.6. We have 

(Jo), = 0, 

which satisfies the third constraint in (|4.13p self- consistently. 
Proof. By definition of the Bloch decomposition, we have 

(5.24) H u n , k = + k f + vo(z)\u n ,k(z) = E n , k u n> k{z). 
Differentiating (|5 . 24|) with respect to fc gives 

(— iV^ + fc)u n ,fc + i?oVfeU n ,fe = ^kE n ,kUn,k + E n ^kUn,k- 

Since Hq is a self-adjoint operator, the above equation taken the inner product 
with u n ^k yields 

(«n,fc|— + fe|tin,fc}ia(r) = Vfc-E„,fc. 

Therefore by ([OH]) . 

(Jo)z = E / ( U ™,fcM V C + fe l u n,(e)L2(r) dfc = E / Vfc£„,fcdfc = 0, 

n<Z n<Z r * 
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where the last equality is due to the periodicity. □ 

Similar to the derivation of (|5 . 1 9[) . by making use of Lemma 15.21 and (|4. 1 3[> . 
direct calculations from (|4.28p give the following proposition. 

Proposition 5.7. The average of Ji(t,x, z) is given by 
(5.25) {J 1 ) z =2Jm J2Y.fi ^^l iV cl u ™^)^(r) e " m " (fc)( ^ T) 
x (^(un.kli^ k\u m ,k) L^(r) ■ V x U (t) + (M„,fe|wi(r)| 

^m,fc/L 2 (r) 

+ (u m ,k\i^7c\ u n,k)L 2 (r) ■ Ao(t)^J dfcdr - A {po)z- 
Then (15.25)) implies ([5T4")) by taking the ^-average of (|4.12|) . 

6. Effective equations in the frequency domain 

We now derive the effective equations in frequency domain. We start with the 
following proposition of the Fourier transform. 

Proposition 6.1. Define the function h(t) = J Q hi(t — r) • h2(r) dr. then 

h = h\H v ■ h2H Vl 
where H v (t) is the Heaviside function oft, 



H v (t) 



1 t > 0, 

otherwise. 



Without loss of generality, we assume that Uo(t) and Ao(t) vanish for t < 0. By 
taking the Fourier transform of (|5.1[) - (|5.4[) and using Proposition 16. 1[ we have 

(6.1) - (S a p + A a p)d Xa d X pUo = B a p{d Xci {A G )p) + p^T tl 

(6.2) - io 2 (A Q ) a - A x (A ) a + (-iu)d Xa U = C a pd X0 Uo + D a0 (A o ) p + J^t- 

The coefficients are given by 



A a p 


— P a j3 ■ 


- Cslw 


xuvy 


f p)zi 


B a f3 


= R a p 




%>vy 


~ 1 9p)z, 


C a /i 


= M a p 


-CV(7- 


- x u vy 




D a /3 


= N a p 


-id* a v(i- 


x-v)- 


_1 9/3)z - <W(po} 2 , 
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where 



M=-X!H / _ 7TT U n,feMm,feK,fe|^l(w)|M m ,fe)l,2 (r) dfe 
n «m>Z' ,r,W + Um,,(fe) 

+ y] T~ 77T<.fc u m.fc(Wn,fe|«i(w)|u m ,fc} i2fr) dfc, 

^) =_ E E T — i 7TT u n,kU* mk (u ntk \iV k \u mik ) L 2 (r) dk 

n<Zm>Z J r* W + W mn (fe) 
+ ~f~ , M <,fcMm,fc(Mn,fcNVfc|^ m , fc ) L2(r) dfc, 

9( W ) = -X1X] / ; , z ; 1 fh \ u ^ku* m , k {u, hk \iv c \u m . k ) L 2 {r) dk 

n< Zm >Z Jr * 

5^ 5^ / 77T<.feWm,fe(Wn,fc|iV C |u m ,fe) L2(r) dfe, 



n<Z m>Z 



and 



= / — ; — - — 7T T (u n!fc |ia fccv |M„ l!fc )(u„. fc |ia fc|3 |u mjfc )dfe 
_ y y / — 7TT( u ™.fe|«9 fea |u„ l!fc }(M„.fc|i9fe 3 |M m ,fc) dfc, 

fia/3(w) = y y / ; 7TT{Un,k\ id k a \ u rn,k){u n ,k\ id C g\ u m,k) dfe 

- y y / — 77T(w„.fc|i9fc Q |w mifc }(u„.fc|ia c |w m .fc)dfc, 

M a/3 (w) = V] V] -f — 7TT(un,k\ id cJ u rn,k)(u n ,k\idkg\u m ,k)dk 

n<Z m>Z 

- y y / 7TT<W ni fe|^ |lt m ,fe)(M„,fe|i9fe |u TO; fe)dfe, 

^H=E E / ; " 7TT( u n,k\ id cJ U rn,k)(Un,k\ id Cg\ U rn,k)dk 

n<z m >z J r- u + u mn (k) 

- y y T — 7Yr(u n ,k\ id C Ju m ,k}(un,k\ id C B \ u m,k)dk. 

We need the following proposition to further simplify the equations. 
Proposition 6.2. 

(r) — tbJmn 
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Proof. Similar to (|5.24|) one has 

H u m . k = + k f + v a {z)^u m . k {z) = E m>k u m>k (z). 

Differentiating it with respect to k gives 

(— + k)u m .k + H Q VkU,n t k — Vfe-B TO) feU m ,fe + E m ,k^kU m ,k- 

Since Hq is a self-adjoint operator, the above equation taken the inner product 
with u Ui k produces 

(U n ,fe|— C,\u m ,k) L 2 {T) = (E m ,k — £ra,fc)(Wn,fc| V 'fe|u m> fc) L 2 {T)i 

which implies the conclusion. □ 
Lemma 6.3. 

(6.3) 23m ^2 E f M «,fc( ;z ) M m,A : ( 2; )( u «^N v fel u m,fe)i 2 (r) dfc = 0, 



(6.4) -23m^ ^ +■ (u n ,fc|zd c Ju m , fe ) 



i 2 (H 

x (« n ,k|i3fe fl |u m: fe} L 2( r ) dfc = {po) z S a f3- 

Proof. Since adding any constant vector to A will not change the system (|3.18[) - 
(|3.20p . the values of p\ and (Ji) z remain the same under the transform Aq — > 
Aq + C v where C v is an arbitrary constant vector. 

Note that we have assumed Aq(£) = for t < 0, then (|5.19|) implies that 

x (u n ,feKVfc|n m ,fe)i2( r ) dfcdr = 0, 
which gives (|6.3[) by making use of Proposition 16.21 



Similarly (|5.25p implies 

-^E e r / r 

' ^ rr J — OO J T * 



e «%n(»)(t-T) 



n<Z m>Z 



x {Mn,fc|j^ o |M m! fc) i2( . r:i (w n! fc|i9fc )3 |u m ^) L 2( r ) dfcdr= (po)z$a/3, 
which produces ()6.4p by making use of Proposition 16.21 □ 
Lemma 6.4. 

g(u) = (-iu)f(uj), N a p(oj) - (p ) z S a i3 = (-iLj)M a p(u), 
M a0 (uj) = -R aP {u), R a p{cu) = (-iu))(P a0 (uj) - P^p), 

where 



T ('u„ ! fc|i9fc Q |u mi fe}(M„,fc|i9fc, 3 |u mi fc)dfc, 

n<Z rn>Z 
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which satisfies P£g — —Pg a . 

Proof. Observe that 

iu> mn (k) . —iu -iuj mn (k) 



w+w mn (fe) w+w mn (k)' w-LJ mn (k) u-uj mn (k)' 

Then it is easy to see that Proposition 16.21 along with (|6.3j) implies g(ui) 
(—iuj)f(u>), and Proposition 16.21 along with f)6.4[) implies N a fj(uj) — (po)z5 a {3 
{-iuj)M a p{u). 

Moreover, Proposition 16.21 also implies 



M a p(uj) = - R a p(u>), 
R a p{u) ={-iu))P a p{(J) 



\idk p \Um,k) dfe 
- f7 ^7 >/ r* 



^ {un,k\idk a \u m , k ){u n ,k\idk ff \u m ,k) dfe 

=H w )(p^( w )-^). 



□ 



Note that 

d Xa d X pU = d X/3 d Xa U , P r a p — —Pp a , 
one knows that the equation (|6.1j) will remain the same if we redefine 

A a0 = P a p - Plp{u) - (?* V(J - X.V)- 1 /^- 
Then Lemma |6 .41 implies 

B a /3 — {~iu)A a p, D a p — (—iw)C a p, C a p — -B a p. 

By defining E = —V x Uo + icuA , B = V x x A , the equations (|6.ip - (|6.2l) along 
with S7 X ■ Aq = produce (|3.29l) - (13.32l) . This completes the derivation of the main 
result in Section [ 



7. Conclusion 

One unsatisfactory aspect of this work is that it is limited to short time scales. 
In fact, the behavior at longer time scales is still very much of a mystery, even from 
the viewpoint of formal asymptotics. This is very unsettling. The main technical 
difficulty is the lack of local charge neutrality and the huge potential generated as 
a result. 

There are other important issues that remain. These include the inclusion of 
spin, the interaction with lattice dynamics, defects, instabilities, etc. 
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